pcor.mat <- function(x,y,z,method="p",na.rm=T){
	x <- c(x)
	y <- c(y)
	
	data <- cbind(x,y,z)
	
	xdata <- data[,c(1,2),drop=FALSE]
	Sxx <- cov(xdata,xdata,m=method)
	
	#xdata <- data[,c(1,2),drop=FALSE]
	zdata <- data[,-c(1,2),drop=FALSE]
	Sxz <- cov(xdata,zdata,m=method)
	
	Szz <- cov(zdata,zdata,m=method)
	#Uses the Moore-Penrose generalized inverse of a matrix.
	Sxx.z <- Sxx - Sxz %*% ginv(Szz) %*% t(Sxz) #C11
	
	rxx.z <- cov2cor(Sxx.z)[1,2]
	rm(data,xdata,zdata,Sxz,Sxx,Sxx.z)
	rxx.z
}
pcor.mat <-cmpfun(pcor.mat)

